/*
 * CTDR_GridField.cpp
 *
 *  Created on: 1 déc. 2010
 *      Author: mbenoit
 */
#include "CTDR_GridField.h"
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
using namespace std;
//ClassImp(CTDR_GridField)

//////////////////////////////////////////////////////////////////////////
// This class represent the CTDR_GridField in the CZT detector. CTDR_GridField 	//
// data are read from a ROOT  file and stored in a TTree.This class pass//
// to the simulation the value of the CTDR_GridField at a position x,y,z,	//
// on demand. This class is similar to the class Detecteur.		//
//////////////////////////////////////////////////////////////////////////


//______________________________________________________________________________
CTDR_GridField::CTDR_GridField(TFile &Fichier, int NdivX, int NdivY, int NdivZ,double LX,double LY, double LZ) :
		CTDR_Field(LX, LY, LZ)
{
// Constructor for the CTDR_GridField
// We retrieve 3 TTree from a file. These data tree contain data of the CTDR_GridField inside a logical pixel of
// our detector. The data are obtained from a finite-emelent analysis software, COMSOL multiphysics, and
// we retreive an array of 100*100*100 for nodes inside the model at constant distance. Each data point
// is assigned a unique ID calculated from its position. It allow us to get the memory adress of a CTDR_GridField value
// by calculating the ID associated to the position where we need the value. We then retreive the CTDR_GridField
// value from the tree using this unique ID pointing to its adress in the tree's memory
// The data trees are created one time , prior to simulation, using Treereader.c ROOT script .
// This allow to modify the script to use CTDR_GridField Data from various sources without touching this code.
// The Storage of the data in .root files allow to get in memory large amount of data, compressed with
// an efficient method that allow a fast acces to data in the tree. S
// parameters: (reference to the ROOT file, number of division in X (100), in Y, in Z, size of the model
// in X,Y,Z)

                   char hash[100];// A character chain to store the hash function
                   NDivx=NdivX;
                   NDivy=NdivY;
                   NDivz=NdivZ;
                   LX=LX;
                   LY=LY;
                   LZ=LZ;
		   //We retreive the TTree from the File
                   TEx=(TNtupleD*)Fichier.Get("Ex");
                   TEy=(TNtupleD*)Fichier.Get("Ey");
                   TEz=(TNtupleD*)Fichier.Get("Ez");
		   //We build the ID->Memory adress Index
                   TEx->BuildIndex("I","0");
                   TEx->SetBranchAddress("x",&x_Ex);
                   TEx->SetBranchAddress("y",&y_Ex);
                   TEx->SetBranchAddress("z",&z_Ex);
                   TEx->SetBranchAddress("E",&E_Ex);
                   TEy->BuildIndex("I","0");
                   TEy->SetBranchAddress("x",&x_Ey);
                   TEy->SetBranchAddress("y",&y_Ey);
                   TEy->SetBranchAddress("z",&z_Ey);
                   TEy->SetBranchAddress("E",&E_Ey);
                   TEz->BuildIndex("I","0");
                   TEz->SetBranchAddress("x",&x_Ez);
                   TEz->SetBranchAddress("y",&y_Ez);
                   TEz->SetBranchAddress("z",&z_Ez);
                   TEz->SetBranchAddress("E",&E_Ez);
		   //We print some info on the tree to be sure it is well loaded
                   TEx->Print();
                   TEy->Print();
                   TEz->Print();
};




//______________________________________________________________________________
void CTDR_GridField::SetChamp(double xx, double yy, double zz)
{
//This function Put in the public variable Ex,Ey,Ez, the CTDR_GridField value for the point xx,yy,zz
    double x, y, z;
     //The model is assumed to be periodic in X,Y
     if((fabs(xx)>LX/2)|(fabs(yy)>LY/2)){
                                   x=fmod(xx+LX/2,LX)-LX/2;
                                   y=fmod(yy+LY/2,LY)-LY/2;
                                   z=zz;
                                   }
     //No CTDR_GridField Data over or under the box
     else if(zz>=LZ){
          x=xx;
          y=yy;
          z=LZ;
          }
     else if(z<0){
          x=xx;
          y=yy;
          z=0;
          }
     else {
          x=xx;
          y=yy;
          z=zz;
          };
     //We interpolate the CTDR_GridField at xx,yy,zz, from the surrounding points in the data tree
     Interpole(x,y,z);
     };


//______________________________________________________________________________
void CTDR_GridField::Interpole(double x, double y, double z)
{
// This function calculate the linear 3D interpolation of the GridField at x,y,z,
//   1. We calculate the ID of the 8 surrounding points to (x,y,z)
//   2. We get the GridField value at he adress corresponding to these ID.
//   3  We compute the interpolation and put the interpolated values in Ex,Ey,Ez
   	;
     //Variable for the CTDR_GridField at the 8 points
     double E000,E010,E100,E110,E001,E011,E101,E111;
     //variable for variable change x,y,z->t,u,v
     double t,u,v;

     /*Ex computation*/
     SetPoints(x,y,z);//We calculate the Id of the 8 surrounding points
     //We retreive the value of the CTDR_GridField at the 8 points using their ID
     Long64_t Adresse=Hash(lowx,lowy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E000=E_Ex;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E010=E_Ex;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E100=E_Ex;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TEx->GetEntryWithIndex(Adresse,0);
     E110=E_Ex;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E001=E_Ex;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E011=E_Ex;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E101=E_Ex;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TEx->GetEntryWithIndex(Adresse,0);
     E111=E_Ex;
     //cout << "V111 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;

     //We do a variable change to locate ourself in a cube around origin
     if(highx!=lowx){t=(x-LX*(lowx/NDivx-0.5))/(LX*(highx/NDivx-0.5)-LX*(lowx/NDivx-0.5));}
     else{t=0;};
     if(highy!=lowy){u=(y-LY*(lowy/NDivy-0.5))/(LY*(highy/NDivy-0.5)-LY*(lowy/NDivy-0.5));}
     else{u=0;};
     if(highz!=lowz){v=(z-LZ*(lowz/NDivz))/(LZ*(highz/NDivz)-LZ*(lowz/NDivz));}
     else{v=0;};

     //We compute the interpolation
     Ex=E000*(1 - t)*(1 - u)*(1 - v) +E100*t*(1 - u)*(1 - v) +
                E010*(1 - t)* u *(1 - v) + E001*(1 - t)*(1 - u)*v +
                E101*t*(1 - u)*v + E011*(1 - t)*u*v + E110*t*u*(1 - v) + E111*t*u*v;

     /*Ey computation*/

     Adresse=Hash(lowx,lowy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E000=E_Ey;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E010=E_Ey;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E100=E_Ey;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TEy->GetEntryWithIndex(Adresse,0);
     E110=E_Ey;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E001=E_Ey;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E011=E_Ey;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E101=E_Ey;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TEy->GetEntryWithIndex(Adresse,0);
     E111=E_Ey;
     //cout << "V111 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     if(highx!=lowx){t=(x-LX*(lowx/NDivx-0.5))/(LX*(highx/NDivx-0.5)-LX*(lowx/NDivx-0.5));}// on calcule une fois seulement
     else{t=0;};
     if(highy!=lowy){u=(y-LY*(lowy/NDivy-0.5))/(LY*(highy/NDivy-0.5)-LY*(lowy/NDivy-0.5));}// on calcule une fois seulementy
     else{u=0;};
     if(highz!=lowz){v=(z-LZ*(lowz/NDivz))/(LZ*(highz/NDivz)-LZ*(lowz/NDivz));}// on calcule une fois seulement
     else{v=0;};
     Ey=E000*(1 - t)*(1 - u)*(1 - v) +E100*t*(1 - u)*(1 - v) +
                E010*(1 - t)* u *(1 - v) + E001*(1 - t)*(1 - u)*v +
                E101*t*(1 - u)*v + E011*(1 - t)*u*v + E110*t*u*(1 - v) + E111*t*u*v;

     /*Ez computation */

     Adresse=Hash(lowx,lowy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E000=E_Ez;
     //cout << "V000 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E010=E_Ez;
     //cout <<"V010 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E100=E_Ez;
     //cout << "V100 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,lowz);
     TEz->GetEntryWithIndex(Adresse,0);
     E110=E_Ez;
     //cout << "V110 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,lowy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E001=E_Ez;
     //cout << "V001 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(lowx,highy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E011=E_Ez;
     //cout <<"V011 "<< x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,lowy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E101=E_Ez;
     //cout << "V101 " << x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     Adresse=Hash(highx,highy,highz);
     TEz->GetEntryWithIndex(Adresse,0);
     E111=E_Ez;
     //cout << "V111 " <<  x_Vx <<' ' << y_Vx << ' ' << z_Vx << ' ' << V_Vx << endl;
     if(highx!=lowx){t=(x-LX*(lowx/NDivx-0.5))/(LX*(highx/NDivx-0.5)-LX*(lowx/NDivx-0.5));}// on calcule une fois seulement
     else{t=0;};
     if(highy!=lowy){u=(y-LY*(lowy/NDivy-0.5))/(LY*(highy/NDivy-0.5)-LY*(lowy/NDivy-0.5));}// on calcule une fois seulementy
     else{u=0;};
     if(highz!=lowz){v=(z-LZ*(lowz/NDivz))/(LZ*(highz/NDivz)-LZ*(lowz/NDivz));}// on calcule une fois seulement
     else{v=0;};
     Ez=E000*(1 - t)*(1 - u)*(1 - v) +E100*t*(1 - u)*(1 - v) +
                E010*(1 - t)* u *(1 - v) + E001*(1 - t)*(1 - u)*v +
                E101*t*(1 - u)*v + E011*(1 - t)*u*v + E110*t*u*(1 - v) + E111*t*u*v;
  ;
};



//______________________________________________________________________________
Long64_t CTDR_GridField::Hash(double x,double y, double z)
{
//This function calculate the ID of a position x,y,z
	Long64_t I;
	I=Long64_t((int)x*(NDivy)*(NDivz) + (int)y*(NDivz) +((int)z));
	return I;
};



//______________________________________________________________________________
void CTDR_GridField::SetPoints(double x , double y ,double z)
{
//This function calculate the index of the 8 surrounding points to (x,y,z)
	lowx=floor((x+LX/2)/(LX/(NDivx-1)));
	highx=ceil((x+LX/2)/(LX/(NDivx-1)));
	lowy=floor((y+LY/2)/(LY/(NDivy-1)));
	highy=ceil((y+LY/2)/(LY/(NDivy-1)));
	lowz=floor((z)/(LZ/(NDivz-1)));
	highz=ceil((z)/(LZ/(NDivz-1)));

	if(lowx>NDivx-1) lowx=NDivx-1;
	if(highx>NDivx-1) highx=NDivx-1;
	if(lowy>NDivy-1) lowy=NDivy-1;
	if(highy>NDivy-1) highy=NDivy-1;
	if(lowz>NDivz-1) lowz=NDivz-1;
	if(highz>NDivz-1) highz=NDivz-1;




     //cout << lowx << " " << lowy << " " << lowz <<" " << highx << " " << highy << " " << highz <<endl;
};



